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ABSTRACT 

This paper is a numerical study of the condensation of the warm neutral medium (WNM) into cold neutral medium (CNM) structures 
under the effect of turbulence and thermal instability. It addresses the specific question of the CNM formation in the physical condition 
of the local interstellar medium (ISM). Using low resolution simulations we explored the impact of the WNM initial density and 
properties of the turbulence (stirring in Fourier with a varying mix of solenoidal and compressive modes) on the cold gas formation 
to identify the parameter space that is compatible with well established observational constraints of the H i in the local ISM. Two set 
of initial conditions which match the observations were selected to produce high resolution simulations (1024 3 ) allowing to study in 
detail the properties of the produced dense structures. 

For typical values of the density, pressure and velocity dispersion of the WNM in the solar neighborhood, the turbulent motions of 
the H i can not provoque the phase transition from WNM to CNM, whatever their amplitude and their distribution in solenoidal and 
compressive modes. On the other hand we show that a quasi-isothermal increase in WNM density of a factor of 2 to 4 is enough to 
induce the phase transition, leading to the transition of about 40 percent of the gas to the cold phase within 1 Myr. Given the observed 
properties of the H i in the local ISM, the WNM and individual CNM structures in the local ISM are sub or transsonic and their 
dynamics are tightly interwoven. The velocity field bears the evidence of subsonic turbulence with a 2D power spectrum following 
the Kolmogorov law as P(k) oc /r 8/3 while the density is highly contrasted with a singificantly shallower power spectrum, reminiscent 
of what is observed in the cold ISM. Supra-thermal line width observed for CNM might be the result of relative velocity between cold 
structures. Finally, the cold structures denser than 5 cm -3 reproduce well the laws M oc L 2 - 25_2 - 28 and cr(|v|) oc 0.5 - 0.8 L 1/3 generally 
observed in molecular clouds. 
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1. Introduction 

Understanding the star formation process remains one of the 
main area of research of modern astrophysics. It is directly 
linked to the way interstellar gas is organized and how dense and 
cold structures form. Turbulence plays a fundamental role here. 
It creates local increases of density that are amplified by gravity 
and that, in some cases, might lead to gravitational contraction. 
At the same time the kinetic energy of turbulent motions needs 
to be dissipated in order for the gravitational collapse to happen. 
The overall star formation process is therefore directly linked to 
the way turbulence shapes the density structure of the ISM and 
how kinetic energy is transferred from large to small scales. 

Molecular clouds, where star formation occurs, exhibit 
supra-thermal line widths and large contr ast of density revealed 
by a log-normal PDF of column density ( Vazquez- Semadeni & 
Gar cia|200ll |Goodman et al.|2009[ |Federrath et al.|2010l|Brunt 
201 0) andby powe r spectra of column density shallower (-2.5, 
Bensch et al.|2001| ) than predicted for subsonic or transonic tur- 
bule nce (-3.66), and than wha t is observed in the diffuse Hi 
( Miville-Desche nes et al.|2003] ). For these reasons they are often 
modeled as isothermal and supersonic turbulent flows. On the 
other hand, because turbulence should dissipate in one dynami- 
cal time if it is not maintained at large scale it is thought that star 



formation might occur within a dynamical (or free fall) time. In 
this scenario the gravitational collapse occurs along the cloud 
formation. Therefore, the initial conditions of the star formation 
process depend on the structure of the cold clouds in the atomic 
phase (Hi). 

Most of the mass of the Galactic interstellar medium (ISM) is 
H i and the re are clear obs ervational evidences that it is thermally 
bi-stable (Dick eyet al.|20 03 ); it consists of two thermally stable 
states at significantly different temperatures (the Cold Neutral 
Medium - CNM - and the Warm Neutral Medium - WNM). This 
is a natural result of the density and temperature dependance of 
the cooling processes at play in the diffuse ISM. Like molecular 
clouds the structure of the H i is highly complex with large den- 
sity contrasts (the CNM fills only about 1 percent of the volume 
but has a density 100 times larger than the WNM) but here the 
combination of tu rbulence and thermal instability (Field 1965; 
Fiel d et aT1 [T969 ) is the main driver of the structure. It is im- 
portant to note that the nature of the density fluctuations in an 
isothermal and supersonic flow is conceptually different than in 
a thermally bi-stable fluid. An isothermal turbulent flow with 
supersonic velocity fluctuations w ill naturally develop strong 
density fluctuations (Ap/p ~ M 2 S , | Vazquez-Semadeni||20 1 2b| . 
In such a flow the density fluctuations are transient, associated 
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with the compression and depression of the turbulent motions. 
In a two-phase medium, once a high density structure is cre- 
ated by a local compression it remains at high density even when 
the gas re-expand because it is in pressure equilibrium with the 
inter-cloud medium. A two-phase medium will therefore pro- 
duce long-lived high density structures contrary to supersonic 
flows in which the fluctuations appear and disappear on dynam- 
ical timescales. Therefore the thermally bistable H i offers a nat- 
ural way of producing highly structured initial conditions for the 
star formation in molecular clouds. 

Because of the complex dynamical processes involved, ded- 
icated numerical simulations are useful to develop insights on 
the exact scenario that leads to the formation of cold struc- 
tures in the Hi. Several studies based on numerical simula- 
tions have explored the thermally bi-sta ble physics of the Hi 
since the early one-dimensional work of |Hennebelle & Peraultj 
(1999 ) who showed the effect of compression in a WNM con- 
verging flow on the triggering of the ph ase transition and the 
production of CNM structures. Similarly Koyama & Inutsuka 



( 2000 2002 ) studied the impact of the propagation of a shock 



wave into WNM. They showed how CNM structures form in the 
post-shock region through the thermal instability which operates 
faster than the duration of the shock (galactic spiral shock or su- 
pernova). In both cases the resulting density contrast between 
the CNM and WNM is ~ 100. To obtain a similar contrast in 
isothermal supersonic turbulence, Mach numbers of the order of 
10 are required. These studies have been extended to two and 
three dimensions simulations of thermally bistable flows with 
turbulence induced by a co nverging flo w ([Audit & Hennebelle 



2005; Hennebelle & Audit 2007; Hennebelle et al.|2007U Audit 
& Hennebellel2010HHeitsch et al.|2005, 2006), by the magne 



torotational instabi lity JPiontek & Ostnk er 2004, [20051(2007} 
Piontek et al.|2009|) or by a driving in Fo urier s pace ( Gazol et al. 
2005[|Gazol & Kim|2010[|Seifried et aLpOTTT ). 

It was established that compression of the WNM in transonic 
converging flows triggers the thermal instability; the increase in 
density can be such that the gas cools rapidly to temperatures of 
a few tens of degrees. The cold gas then fragments into clumps, 
in particular through the Kelvin-Helmholtz instability (Heitsch 
|et al.|2006| ). The thermal pressure of the cold gas is in approxi- 
mate equilibrium with the sum of the ram and thermal pressure 
of the surrounding warm gas. Once pushed at high enough densi- 
ties, the cold structures can have densities and temperature close 



to values typical of molecu lar clouds ( Vazquez-S emadeni et al. 
2006||Banerjee et al.|2009) 



In general previous works showed that the efficiency of the 
formation of the CNM depends directly on the nature and ampli- 
tude of the injected energy. One question we wanted to explore 
in this study is related to the driving of the compression. Most 
studies have explored the case of converging flows where the 
increase of the WNM density needed to trigger the thermal in- 
stability is obtained by a sustained ram pressure. The origin of 
such large scale convergent flows is unclear so here we want to 
see if compression in turbulent motions alone, driven at large 
scales, could lead to the production of cold gas. In more details, 
we wanted to see if turbulent motions, at the level they are ob- 
served in the WNM, can induce the phase transition and if not, 
what are the required physical conditions of the WNM to trig- 



ger the formation of cold structures. To do so, like Gazol et al. 
(2005); |Gazol & Kirn] J20l0) ; |Seifried et all <2011) , we used 



turbulen t forcing in Fourier sp ace. In addition, similarly to the 
study of |Federrath et al. ( 2010| ) based on isothermal simulations 
of molecular clouds, we wanted to study the impact of compres- 
sive and solenoidal modes in the injected turbulent field. 



In this paper we present a parameter study of WNM initial 
density and turbulent forcing properties on the formation of the 
CNM in forced turbulence and thermally bi-stable flows. In a 
similar study, [Seifried et aL] ( [201 1| ) chose initial densities (1 to 
3 cm -3 ) that are in the thermally unstable regime, to expect the 
production of a two-phase medium. Our intent is to expand the 
parameter space to identify the dynamical conditions that lead 
to the formation of the cold phase. In addition, one fundamental 
aspect of the current study is that we guide our parameter search 
based on all available observational constraints of the local ISM. 
We present 90 simulations performed on a 128 3 grid allowing 
to identify more precisely the physical conditions that favor the 
formation of cold structures while reproducing in detail the ob- 
served properties of the H i in the local ISM. For two specific 
cases that match the observations we performed two 1024 3 sim- 
ulations in order to study in greater details the properties of the 
H i: temperature, pressure and density histograms, power spec- 
tra, power laws of the cold structures. 

The paper is organized as follow. HERACLES, the numeri- 
cal code used in this study is briefly presented in § [2] The ob- 
servational constrains used to lead the parameter study are de- 
tailed in § [3] The methodology is presented in § [4] The results of 
the parametric study and of the high resolution simulations are 
presented in § [5] and § [6j and then discussed in § [7] The main 
conclusions of the study are given in § [8] 

2. Numerical methods 

Despite their great importance in the interstellar medium, we 
won't consider magnetic field or gravity in this study to focus on 
the impact of vortical and compressive turbulence motions on the 
formation of CNM. All the following simulations are thus hy- 
drodynamical. We performed our simulations with HERACLES 
( [Gonzalez et al.||2007] ). The numerical methods are similar to 
those described in paper from | Audit & Hennebelle (2005 ). 



2.1. Euler equations and Godunov scheme 

The Euler equations for a radiatively cooling gas are the clas- 
sical equations of hydrodynamics in order to treat the flows at 
high Reynolds number, which is the case in the ISM. An exter- 
nal force f is adding to generate the turbulent motions. It has the 
dimension of an acceleration. We also include cooling and heat- 
ings terms based on the prescription of Wolfire et al. (2003} (see 
Fig.[T]) in the net cooling rate X. 



dtP + V • [pu] = 0, 

d t [pu] + V • [pu <S> u + P] = pf , 

d t E + V • [u(E + P)] = pi • u - £(p, T) 



(1) 

(2) 
(3) 



where p is the mass density, u the velocity, P the pressure and E 
the total energy. The considered atomic gas is perfect and ideal 
with an adiabatic index y = 5/3 and a mean atomic weight 
ji = I Amu (wr is the mass of the proton). The weight \i takes 
into account the cosmic abundance of other elements : helium 
and metals, and so the mass density p = /uxn. The net cooling 
function is defined as L - n 2 A - nT (expressed in erg cm -3 s" 1 ), 
where A is the cooling and T the heating. They include the most 
relevant processes in the diffuse medium for our range of tem- 
perature [10K, 1000K] and are plotted and quickly described on 
the Figure [T] (for more details, see | Audit~& Hennebelle ]|2005 



The equations solver in HERACLES is a second order Godunov 
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Fig. 1. Heating and cooling proce s ses 
HERACLES based on IWolfire et all 420031. 



implemented 
We 



represent 

the respective energy densities versus the temperature at the 
fixed density n=l cm -3 , typical of the WNM. The total cooling 
(black solid curve) is decomposed in its different components: 
CII (dark blue), OI (light blue), recombination on interstellar 
grains (red), Lyman a (green). The dashed blue line is the 
heating dominated by the photo-electric effect on small dust 
grains. We considered a spatially uniform radiation field with 
the spectrum and inte nsity of the Habing field (G o/ 1.7 where 
G is the Draine flux) ( [Habing 1 1 968} |Draine| 1 978] ) . 



scheme (see Audit & Hennebelle 2010) for the conservative part. 
It is a MHD solver of type MUSC L-Handcook( |Londrillo & Del 
Zanna 2000; Fromang et al.|2006| ). The cooling is applied after 
(see |Audit & Hennebelle|2005| ). 

2.2. Turbulence forcing 

The force field f generates the large scales motions and is applied 
in Fourier space. We define the amplitude given to the field as 
vs . The magnitude of the applied force is thus F s = v 2 s /L s where 
Ls = LboxAo is the characteristic integration scale and ko = 2 the 
characteristic wave number of the stirring. The auto-correlation 
time of the process is given by T$ = (L s /Fs) 1/2 = ^s/^s- The 
stirring is applied at large scale, and thus to the small wave 
numbers k that fulfill < \k\ < 2& . The modes of the turbulent 
forcing are computed with the stochastic process of Qr nstein- 
Uhlenbeck ( [Eswaran & Pope|| 19881 ISchmidt et al.|[2006]). The 
method is similar to |Schmidt et al.[ ( |2009| ) and |Federrath et al.| 
2010| ). A Helmholtz decomposition is applied to the random 



field and project it on its compressive and solenoidal compo- 
nents: 



01% 



&u + (i 



2f # 



(4) 



Choosing a spectral weight £ = 1 creates a physical force di- 
vergence free, i.e. purely solenoidal. Taking £ < 1 generates di- 
latational components which compress or rarify the gas. Finally, 
f = creates a rotation-free field, i.e. purely compressive. In a 
natural state, meaning that the energy is naturally distributed be- 
tween the compressive and solenoidal modes as the 2:1 mixing, 
f takes the value of 0.5. In this case, the Helmholtz operator is 
proportional to the identity operator and the projection does not 
alter the energy repartition between the different modes. 



3. Observational constraints 

In this section we present global properties of the Galactic H i 
as deduced from the observations. These properties guided us in 
the choice of the initial conditions of the simulations and on the 
qualification of the results obtained. The simulations presented 
in this study are dedicated to the H i thermal instability and to the 
formation of CNM out of pure WNM gas. What are reasonable 
properties of the WNM and what are exactly the observational 
constraints ? 



3.1. Pressure and thermal bi-stability of the H i 

Because of the density and temperature dependance of the heat- 
ing and cooling processe s in the diffuse ISM (see Fig. [T] ) it has 
been shown early on b y |Field| (| 1 965| ) ; [FTeld et al.|p96 9) and up- 
dated by Wolfir eet al.| (p~995 , 2003 ) that there is a range of pres- 
sure where the neutral atomic gas can be in two thermally stable 
phases in pressure equilibrium. This is generally displayed in di- 
agrams like the ones shown in Fig. [2] where the solid curves in 
the two panels indicate the locus of thermal equilibrium in the 
P — n and T - n parameter spaces. 

The thermal stability is not preserved everywhere along these 
curve but only where dP/dn > 0. The two dash-dotted verti- 
cal lines mark the boundaries in density where the Hi is ther- 
mally unstable. Any perturbation of the gas sitting on the equilib- 
rium curve in the thermally unstable range would move towards 
the cold or warm branches. From Fig. [2] one can appreciate the 
steepness of the equilibrium curve in T - n space, highlighting 
the great difference in temperature and density of the two stable 
phases of the H i. 

It is important to point out that the bi-stability of the Hi 
is only present in a given range of pressure, determined by 
the maximum density of the warm branch (n 
the minimum density of the cold branch in 
P < 1000 K cm" 3 there is no CNM and for P 
all the H i is cold. 

On the top panel the dashed curve gives the average pressure 



< 0.8 cm" 3 ) and 
- 7 cm" 3 ). For 
> 6000 K cm" 3 



observed in the local ISM. According to Jenkins & Tripp (201 1 



the pressure of the cold phase follows a log-normal distribu- 
tion with an average of log 10 (P/&) equals to 3.58 and minimum 
and maximum values (at the 99 percent level) equal to 3.04 and 
4.1 1. The minimum CNM pressure reported by Jenkins & Tripp 
( 201 1) corresponds exactly to the minimum CNM pressure (3.06 
dex) allowed by the thermal instability (see Fig. [2]) as predicted 
by the mode ling of the heating and cooling processes of Wolfire 
et al. ( 2003|). On the other hand, the maximum pressure reported 
by | Jenkins & Tr ipp ( 201 1 ) is significantly higher than the maxi- 
mum pressure allowed for the thermally stable WNM (3.74). In 
fact, the mass fraction of the CNM found at higher pressure by 
Jenkins & Tripp ( 2011) is about 30 percent. In those cases most 
of the WNM is in the thermally unstable regime or in a cooling 
phase; in any case the warm gas is an over-pressured state that is 
likely to be transient. 



3.2. WNM density and temperature 

Here we want to estimate what is the true density (as opposed 
to space-averaged) and temperature of the WNM in the local 
ISM. Typical numbers cited are n ~ 0.2 - 0.5 cm" 3 and T ~ 
6 000 - 10 4 K (|Ferriere||200lt , but, as emphasized by |Kalberla 
Pede s (2008), the density of the WNM varies greatly with 
galacto-centric radius and galactic height z. It is also subject to 
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Fig. 2. Equilibrium temperature and pressure as a function of 
density. These equilibrium curves were c omputed equating th e 
heating and cooling processes described in Wol fire et al.| (2003). 
The dashed line correspon ds to the average pressur e of the cold 
diffuse ISM measured by Jenkins & Tripp| ( |2011| ). The dotted 
lines show the maximum and minimum pressures that include 99 
percent of the data of Jenk ins & Tripp| ( [2011| ). The dash-dotted 
lines show the density range of the thermally unstable phase. 



Fig. 3. Vertical distribution of the WNM density, temperature 
and pressure. The density profile assum es n(z = 0) = 0.5 cm" 3 , 
in accordance with the mean pressure of Jenkins & Tripp ( 201 1) 
and th e heating and cooling processes described in Wolfire et al. 
( 200 3). The scal e height of the WNM is the value of |Dickey & 



Lockman (1990) : HWHM = 265 pc. The temperature is com- 



puted at each position z assuming thermal equilibrium at density 
n(z). The pressure is simply P/k = nT. 



1.00 



local variations d ue to star formation activity (e.g. supernovae, 
Joung et al.|2009| ) and between the arm and inter- arm regions. 
According topickey & Lockman (1990) (see also 



Celnik 



|eTaLl[T979l |Malhotra||19951|Kalberla et al.||2007| ), the distribu- 
tion of the H i density with z close to the Sun can be decomposed 
in a narrow Gaussian of HWHM=106pc, a larger Gaussian 
of HWHM=265 pc and an exponential with a scale height of 
403 pc. Generally the narrow Gaussian is attributed to the CNM 
and the sum of the two others to the WNM. In this scheme, and 



still following Dic key & Lo ckman 



1990| ), the space-average 
density of the CNM and WNM at z = are 0.395 cm -3 and 
0.172 cm" 3 respectively. These values are weighted by the vol- 
ume filling factor / and thus correspond to n f, where n is the 
true volume density. The filling factor of the WNM is difficult to 
estimate; according to Kalberla & Kerp (2009) 



it is ~ 0.4 which 



implies a true density of ^wnm = 0.43 cm J . 

This value is close to the one that can be derived based on 
the knowledge of the heating and cooling proc esses. For the av 
erage pressure of P = 3980 K cm -3 found by ( [Jenkins & Tripp 
|2011|), using the heating and cooling curves of Wolfire et a] 



( 2003) and assuming a spatially uniform radiation field with the 
spectrum and inte nsity of the Habing field (Go / 1.7 where Go is 
the Draine flux - Habing 1968; Draine 1978), the equilibrium 
density and temperature of the WNM are n wnm = 0.5 cm -3 and 
r wnm = 7960K. 

Fig. [3] shows how n, T and P vary with z assuming a 
scale height of HWHM=265 pc for the WNM disk at R = R Q 



( [Dickey & Lockman|1990HMalhotra|1995l|Kalberla et al.|2007) . 
Assuming a true volume density tzwnm = 0.5 cm" j at z = 0, the 
corresponding temperature profile T(z) of the WNM can be esti- 
mated using the thermal equilibrium cur ve based on th e heating 
and cooling processes described in Wolfire et al. (2003). It shows 
that the variations with z of the WNM properties due to the hy- 
drostatic equilibrium become important only at scales larger than 
100 pc. 
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Fig. 5. Temperature histogram weighted by density (on top) and 
density pdf (at the bottom) for simulations with the following 
initial conditions: n = 1 cm -3 , £=0.2, vs = 12.5 km s" 1 , and dif- 
ferent resolutions: 512 3 (solid black line), 256 3 (dotted blue line) 
and 128 3 cells (dashed green line) at 40 Myr. 
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time (Myear) time (Myear) 



Fig. 4. Time evolution of the CNM mass fraction /cnm, the mean velocity dispersion cr turb and the mean Mach number At bs for 
simulations with the following initial conditions: on the left no = 1cm" 3 , £=0.2, vs = 12.5 km s" 1 , on the right n=2 cm" 3 , £=0.5, 
vs = 7.5 km s" 1 , and different resolutions: 512 3 (solid black line), 256 3 (dotted blue line) and 128 3 cells (dashed green line). 



3.3. WNM turbulent velocity dispersion and Mach number 

Apart from its density and temperature, the other relevant pa- 
rameter for this study is the level of turbulence in the WNM, or 
equivalently its Mach number. The direct estimate of the Mach 
number of the WNM is difficult due to the fact that its temper- 
ature can not be easily measured. The traditional method to es- 
timate the Hi temperature is based on a comparison of 21 cm 
absorption and emission measurements. Due to its low opac- 
ity, the warm phase produces almost undetectable absorption 
once observed in front of a strong radio-continuum source. Only 
a few detections were reported so far (e.g. Carilli et"aL||1998| 
Kan ekar et al.|2003||Begum et al.|2010"| ), with deduced tempera- 
tures closer to the thermally unstable regime (T ~ 500 - 6000 K) 
than to what is expected for the WNM (T ~ 8000 K). 

Nevertheless there are indirect evidences that the warm 
ISM is subsonic. Recent findings on the warm ionized medium 
(WIM) by Gaensler et al. ( 2011|) confirmed the wor k of 
Armstrong et al.| ( |1995"] Tand |Redfield & Linskyj p004| ) that 
showed that the WIM has the properties of a low Mach num- 
ber and incompressible turbulent flow. Regarding the WNM, an 
indication of its low Mach number comes from the smoothness 
of the 21 cm profiles. |Haud & Kalberla| ( [2007]) found that the 
simplest 21 cm spectra of the LAB data ( [Kalberla et al.||2005] ) 
at high Galactic latitudes are well modeled by a single Gaussian 
of <x tot = 10.2 ± 0.3 km s" 1 . Here we assume this value to be 
representative of the WNM total velocity dispersion in the Solar 
neighborhood. 

The WNM velocity dispersion is the quadratic sum of the 
thermal (<x t h e rm) and turbulent (cr turb ) contribution to the gas mo- 
tions, integrated along the line of sight. Considering the specific 



line of sight at the Galactic pole, the thermal contribution to the 
line width of the WNM is an integral over Galactic height z' 



k f L n(z)T(z)dz 



^ n(z)dz 



1/2 



(5) 



where n(z) is assumed to be a Gaussian with HWHM=265 pc 

( Dickey & Lockman 1990 ). Using the corresponding T(z) shown 
^ c; ^ nl ;„f^„^t;^,v c Pel ^ _ o „-! « , ™ 



in Fig[3J integrating Eq.|5| gives <x t herm = 8.3 km s~ , a value es- 
sentially identical to the thermal width of an isothermal 8000 K 
gas (cr therm = 8.1 km s" 1 ) which reflects the fact that the WNM 
temperature does not vary signi ficantly with Galactic height. 
Assuming cr tot = 10.2 km s" 1 (|Haud & Kalberla||2007]), the 



contribution of turbulent motions to the observed line width is 
0"turb = 5.9 km s -1 . The fact that cr tmb < cr therm is another indi- 
cation that thermal motions dominates over turbulent motions at 
any scale in the Solar neighborhood. 

Due to the nature of the turbulent cascade, the turbulent ve- 
locity dispersion increases with scale like <x tur b = o~(l)l q , where 
ci ~ 1/3 for subsoni c turbulence and, following the notation of 
Wolfire et al. (2003), <x(l) is the velocity dispersion (in kms" 1 ) 
at a scale of 1 pc and / is the scale (in pc). Assuming that the tur- 
bulent contribution to the observed line width is representative 
of motions at the scale of the H i scale height at the Sun location 
(I = h z = 265 pc), one obtains <x(l) = 0.92 km s" 1 . This calcula- 
tion is similar to the one of [Wolfire et al.| (|2003) who estimated 
<t(1) = 1.4 km s" 1 for the WNM. The fact that we obtain a sig- 
nificantly lower value can be attributed mostly to the fact that we 
took into account the thermal contribution to the line width. 
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3.4. CNM mass fraction and volume filling factor 

One important quantity we want to reproduce in the simulations 
is the amount of mass in the cold phase. The CNM mass frac- 

(2003) 



tion was estimated to 40 percent by |Heiles & Troland 



by comparing 21 cm absorption and emission measurements on 
79 random lines-of- sights. This is also compatible w ith the de- 
composition of the H i density distribution with z of Dickey & 
Lockman|( |1990| ); their narrow component (h - 106 pc) that can 
be associated to the CNM contributes to 44 percent of the total 
H i surface density of the disk. Because the scale height of the 
CNM is smaller than that of the WNM, the CNM mass fraction 
varies with z. The CNM contributes to 60 percent of the H i mass 
forz< 200pc. 

Like for the WNM one can compute the expected CNM vol- 
ume dens ity based on the heating and cooling processes ( Wolfire 
et al. |2003] ). For the average pressure of Jenkins & Tripp ( 201 1 ) 
(P = 3980 K cm -3 ) the equilibrium density and temperatur e are 



cnm 



= 111 cm" 3 , T„ 



35.9 K. Following Dickey & Lockman 



(1990), the space-average density of the CNM is ~ 0.4 cm" J at 
z = 0, which implies a volume filling factor of the cold phase 
lower than 1%. 



4. Methodology 

Based on the observational constraints described in the previ- 
ous section, we conclude that the WNM in the Solar neighbor- 
hood has a mean temperature T ~ 8000 K, a mean density of 
n ~ 0.2 - 0.5 cm -3 and subsonic turbulent motions following 
0"turb = 0.92 (//lpc) 1/3 (kms -1 ). The first step of this study is 
to evaluate what are the plausible physical conditions that lead 
to the formation of CNM gas, of the order of 40 percent in 
mass, from WNM gas with such physical properties. We ex- 
plored those conditions with 90 simulations at low resolution 
(128 3 cells). In this section we validate this choice by testing 
the use of 128 3 simulations for the exploration of the parameter 
space by comparing with higher resolution (256 3 and 512 3 ) sim- 
ulations, after having described the characteristic scales involved 
in the simulations. Two sets of initial conditions that match the 
observations will then be used to produce high resolution simu- 
lations (1024 3 ). 



4.1. Characteristic scales 

4.1 .1 . Dynamical time and cooling scale of the WNM 

If, during the perturbation of an initial WNM cloud, the cooling 
processes act slower than the compressive ones, the gas will tend 
to return to its initial state of temperature and density by decom- 
pression. The transition WNM-CNM becomes possible when 
the radiative time is lower than the dynamical time (Hennebelle 
& Perault|1999| ), allowing the condensed gas to keep its higher 
density even during the decompression: 

(6) 



^cool ^ ^dyn- 

The dynamical time is defined as the time needed by an atom to 
cross the cloud (of size d) at the sound speed velocity c$\ 

d 

%n - — • (7) 

The radiative time is the characteristic time needed by the 
medium to cool. It is thus equal to the ratio of the energy (U) 
and the energy transfer (e) applied to the medium: 

_ U 

^cool = — • (8) 
6 



We thus can calculate the critical size below which the WNM 
will be unstable and will condense into the unstable regime or 
CNM structures. We consider typical conditions of the WNM 
(T = 8000 K and n - 0.5 cm -3 ) and a static medium, assuming 
that the energy is fully thermal. Using equation [6] and the heat 
capacity per mass unity Cy = k B /mn(y - 1): 



U = muCyT the gas energy and 
e = nA the cooling applied to the WNM 

The cooling time thus becomes: 



^cool — 



k B T 



(9) 
(10) 



(11) 



(y - \)nK 

Considering the main cooling processes of the diffuse Hi de- 
to evaluate A and cs = ^jyk^T 7(1.4 mn), 



3.1 



scribed in Section 

one obtains for the WNM cooling scale: 

^cool < tdyn => d > /i coo l = ^cool x <?S,WNM ~ 22 pc. (12) 

Thus, for the condition t coo \ < t^ yn to be observed and the ther- 
mal instability to develop, the initial scale d of the perturbed 
zone must be greater than A coo \ ~ 22 pc. Similarly, a perturba- 
tion of size lower than A coo \ do not efficiently trigger the phase 
transition. 



4.1 .2. Sizes of the simulations 

To study the WNM-CNM transition it is important to evaluate 
the scales of cold structures to adjust the numerical resolution. 
To first order one can estimate the size of a cold structure result- 
ing from the condensation of a WNM cloud with a typical size 
of 22 pc. Because the CNM density is a hundred times higher 
than the WNM density, the CNM structures should be a hun- 
dred times smaller than the initial WNM cloud, leading to a 
rough size of 0.2 pc. In fact, as illustrated by Hennebelle & Audit 
( 2007]), the size of the CNM clumps follow a distribution with 
a maximal size of 0.6 pc. Moreover, the Field length A? (Field 
1965; Begelman & McKee 1990), which give the typical scale 
of the thermal front between the WNM and the CNM, and the 
front shock in the CNM are of the order of 10" 3 pc. In principle, 
simulations of the thermally bistable Hi require a scale range 



from 10 3 to ~ 20 pc but, as suggested by Hennebelle & Audit 



(2007), the important scale to resolve in simulations of thermally 
bistable and turbulent Hi is the 'sound crossing scale' which is 
the sound speed times the cooling time. As mentioned earlier, for 
typical values of the density, this scale is ~ 22 pc and ~ 0.2 pc in 
the WNM and CNM respectively. 

We chose to assure the phase transition with a box size equal 
to 40 pc and therefore larger than A coo \. In this study we used first 
a box size of 128 3 , leading to a resolution of 0.3 pc. Even though 
we barely resolve the typical cold structure scale, using higher 
resolution simulations we will demonstrate that this set up can 
be used to explore efficiently the parameter space for the for- 
mation of the CNM. The second part of the results presented in 
this paper is dedicated to the detailed study of 1024 3 simulations 
(resolution of 0.04 pc) for specific initial conditions. 

4.2. Effect of resolution 

Here we validate the use of low resolution simulations in the 
context of a parameter exploration study. Figure [4] shows the 
evolution over time of the CNM massive fraction, the disper- 
sion of the line of sight turbulent velocity <x tur b and the ob- 
servational Mach number for two representative simulations 



6 



Saury et al.: The structure of the thermally bistable and turbulent atomic gas in the local interstellar medium 






10 . 




8 - 


C/) 




E 


6 L 






§ 






4 - 






'co 






2 - 




0" 








5 km/s 
7.5 km/s 
1 km/s 
12.5km/s 
1 5 km/s 
20 km/s 




4 6 

n (cm 3 ) 



10 



Fig. 6. Behavior of /cnm ^£ and <x tur b (on the right) with the variation of the initial condition no, for the different values of 
the large scale velocity vs (from 5 to 20 km s" 1 ). The dotted lines frame the observational values of /cnm and <T tur b that we want to 
reproduce. 
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Fig. 7. Maps of the logarithm of the column density (density integrated along the direction z). The initial density increases from the 
left to the right and takes the following values: 1.0, 1.5, 3.0, 10.0 cm" 3 . The projection weight and the large scale velocity are fixed 
to £ = 0.5 et vs = 7.5 km s" 1 . All maps are displayed with the same column density range [0.6 10 20 cm -2 , 40 10 20 cm -2 ]. 



([n=l cm" 3 , £=0.2, v s = 12.5 km s" 1 ] and [n=2cm" 3 , £=0.5, 
v s = 7.5 km s" 1 ]) done on 128 3 , 256 3 and 512 3 grids. 

We calculate the cold mass fraction /cnm as the sum of the 
density in the pixels where T < 200 K over the total density in 
the box. The turbulent velocity dispersion is computed the fol- 
lowing way (see |Miville^D eschenes & Martin 2007 , equation 8): 



9 2 



Lbox 

z 



9 x n(x,y,z')v 2 z (x,y,z') 



Lbo x n(x,y,z>) 



-C 2 (x,y) 



where C(x, y) is the velocity centroid map: 

C(x,y) 



Zzt°o X n(x,y,z')v z (x,y,z') 



Z?:o X n(x,y,z>) 



(13) 



(14) 



In numerical simulations, the Mach number is often com- 
puted as the ratio of the velocity and the sound speed on each 
cell: 



Mheo =< |U|/C S >. 



x,y,z 



(15) 



where cs = yj(yP(x, y, z))/p(x, y, z). Because this parameter can 
not be directly deduced from observations we consider an ex- 
pression for the Mach number calculated from integrated quan- 
tities along the line of sight: the turbulent velocity dispersion and 
the thermal velocity dispersion: 



M 



obs 



0"turb 

= < >x 

^"therm 



(16) 



where the thermal velocity dispersion is defined as 



h ^n{x,y,z')T(x,y,z') 



therm 



Obviously the totaQ velocity dispersion is simply 



2 _ 2 2 
^"tot _ °"turb + ^"therm- 



(17) 



(18) 



All quantities shown in Fig [4] are spaced averaged over the 
box at different time steps. The mean values of cr turb and At bs 
are in the same range, and their variations are of the same order, 
as illustrated in Figure [4] These two quantities are thus statisti- 
cally well represented even in small resolution simulations. 

The temperature and density pdfs (Fig. [5]) and time evolution 
(Fig. [4]) show that the ratio of cold gas depends on the numer- 
ical resolution. Higher resolution runs lead to a slightly higher 
amount of cold gas. The cold peak of the temperature histogram 
(Fig. [5]) is lower in the 128 3 and 25 6 3 cases. The time evolution 
of /cnm is increasing slower for the smallest resolution run and 
/cnm is then lower. However the massive fractions of gas with 
T < 200K reached at 40 Myr are similar (0.27, 0.31 and 0.33). 
We then deduce that the estimation of /cnm on 128 3 simulations 
is still worthwhile if we remember that it increases slightly with 
resolution. 

The same resolution study on simulations with different 
initial conditions presents similar results. Like Seifried et al. 



or observed in the case of an optically thin line 
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Fig. 8. Velocity dispersions histograms: on top cr therm , in the 
middle <x tur b and at the bottom cr tot , for different temperature 
ranges: in black the all box, in dark blue the WNM (T > 
5000 K), in blue the unstable gas (200 K < T < 5000 K) and 
in light blue the CNM (T < 200 K) for the following high reso- 
lution simulation: 1024 3 , no=2.0 cm -3 , f = 0.5 et vs=7.5 km s" 1 . 



( |2011| ), we conclude that the resolution does not affect the time 
evolution of <x tur b and At bs and has a moderate effect on the 
formation of CNM. Therefore doing a statistical study on small 
simulations (128 3 ) is legitimate. Higher resolution simulations 
are needed to study in details the CNM properties, especially its 
structure. 
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Fig. 9. Density, temperature and velocity (z-component) pro- 
files for a given simulation (1024 3 , no=2.0 cm -3 , £ = 0.5 et 
vs=7.5 kms" 1 ) and for two distinct lines of sight characterized 
by their value of cr turb (CNM). On top: <x turb (CNM) =0.14 kms" 1 
and at the bottom: <T tur b(CNM)=4.7 kms -1 . 



4.3. Convergence 

In order to evaluate the time needed for the 128 3 simulations 
to reach a stationary state, we evaluated cr tmh , At bs and /cnm 
at each Myr. Most simulations converge in 10 Myr but some 
needed a much longer time to transit to cold phase and never 
achieve a stationary state for the massive fraction of CNM. 
However we choose to stop all simulations at 40 Myr because of 
the typical time between two su pernovae explosions estimated 
to be around 20 Myr (see Ferriere|2001 equation 17 at z = in 
a volume of (40 pc) 3 ). 



5. Parametric study 

According to Fig. [3] we can consider constant n and T as ade- 
quate initial conditions for a box size of 40 pc centered on z = 0. 
For all simulations the initial conditions are a uniform and static 
WNM at a temperature of T = 8000 K. The boundary conditions 
are periodic. The velocity field is modified at every time step by 
the turbulent stirring in Fourier space. 

This section presents the parameter study we did on the ini- 
tial density and the two properties of the stirring velocity field: 
the spectral weight £ and the large scale velocity vs. We per- 
formed 90 simulations with 128 3 cells, varying the density from 
its fiducial WNM value 0.2 to 10 cm -3 , the velocity at large scale 
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vs from 5 to 20 km s" 1 and the spectral weight from (rotation 
free field) to 1 (divergence free). We computed for each simula- 
tion the CNM mass fraction /cnm, the mean Mach number At bs 
(Eq. 16) and the mean velocity dispersion (Eq. 13 > at each Myr 
and averaged them over time after a stationary state has been 
reached, that is on the time interval [20 Myr,40 Myr]. 

5.1. Effect of the initial density increase 

Tableland Figure|6]present the results for the 48 simulations for 
which the turbulence has been kept to a natural state, meaning 
that the spectral weight is fixed to f = 0.5 and that the Helmholtz 
projection operator (Eq|4]) is equal to half the identity and there- 
fore keeps the rate of compressible over solenoidal modes of the 
velocity field in its initial state. In this case, we only vary the 
large scale velocity vs and the density. Obviously, as the initial 
temperature is 8000 K for each simulation, increasing the initial 
density is equivalent to increasing the initial pressure. 

We first observe that simulations with initial values of the 
WNM density between 0.2 and 1.0 cm -3 do not transit to cold 
gas. At larger densities, the efficiency of the CNM formation in- 
creases rapidly with the density. For no = 0.2 or 1.0 cm" 3 , the 
cold mass fraction is indeed very close to zero while a small in- 
crease of the density (no = 1.5 cm" 3 ) leads to /cnm= 0-3 - 0.4 
and densities greater then 3.0 cm" 3 allow more than half of the 
mass to transit (see Fig(6] left part). This behavior is also well 
represented by the differences regarding the structures on the 
integrated density maps (Figure|7|). The map on the left, that 
only contains WNM, is completely smooth while the structures 
appear clearly when the initial density increases. One can no- 
tice that the equivalent pressure at no = 1.5 cm" 3 is equal to 
12000 K cm" 3 , much larger than the maximum pressure allowed 
for the WNM (see Sec. |3.1| ). The gas starts in the thermally un- 
stable regime; at these values of density and pressure the gas is 
indeed located above the stable branch of the WNM and thus 
in a cooling zone that allows it to evolve rapidly and efficiently 
towards the cold gas. 

From Table [T] and Figure|6]one could note that the final tur- 
bulent velocity dispersion is lower for higher initial densities; 
cr turb is thus lower for higher values of the cold mass fraction 
/cnm- To understand this behavior we studied the distributions 
of the three following velocity dispersions: cr turb , <x t herm and cr tot 
computed for each thermal phase in a given simulation at high 
resolution (1024 3 ,^ = 2.0cm" 3 , f = 0.5 and v s =7.5 kms" 1 ). 
All distributions are plotted on Figure [8] As expected <x t herm de- 
pends on the temperature of the gas. We note that it is also the 
case for the distributions of cr turb , the one computed in the CNM 
phase (T < 200 K) being especially different from the other two 
thermal phases. Its peak lies at very small values between 0.1 
and 0.2 kms" 1 and presents a wide tail reaching 5 kms" 1 . To 
better understand the shape of this distribution, we looked at 
two distinctive lines of sight for which <T tur b(CNM) is equal to 
0.14 and 4.7 km respectively. The density, temperature and ve- 
locity (z-component) profiles are plotted on Figure [9l The line 
of sight with a low velocity dispersion (top of Fig]9|) is com- 
posed of four CNM structures with temperatures between 100 
and 500 K. These clumps have homogeneous internal velocities 
and their relative velocity is small (~ 3km s" 1 ), explaining the 
low value of <x tur b(CNM) integrated along the line of sight. On 
the other hand, if the internal velocities of the CNM structures 
on the second line of sight are also homogeneous, their relative 
velocity is high and close to 15 kms" 1 . Thus, the shape of the 
histograms of cr turb (CNM) reveals that individual clumps have an 
internal velocity dispersion close to 0.2 km s" 1 and that the large 



values are the result of the relative velocities between clumps. 
Besides, the three graphs on Figure [8] show that the shape of the 
<x tot histograms is dominated by th e thermal motions sugges ting 
a subsonic Mach number. Like jHeitsch et al.| ( |2005| |2006| ) we 
conclude that the cold gas structures have subsonic internal mo- 
tions and tha t the supersonic Mach nu mber generally observed 
in the CNM ( [Heiles & Troland||2003| ) can be the result of their 
relative motions. 

A similar explanation is proposed for the increase of the ob- 
served Mach number with density (see Table [TJ. By increasing 
the density, the cold mass fraction increases and so the number 
of cold clumps in the box. These have relative velocities inher- 
ited from the WNM and therefore supersonic with regard to their 
low temperature. 

As mentioned by Audit & Hennebelle ( 2005), we note that 
at a given initial density, the efficiency of the CNM formation 
decreases when the large scale turbulent velocity increases (see 
Fig(6| on the left) and thus when the Mach number increases. 
Vazquez-Semadeni (2012a) argues that supersonic turbulence 
acts faster than TI and thus dominates the evolution of the cloud. 
It has been recently confirmed by Wal ch et aL](|2Ql lj) . As st ated 
in Eq.[6]and as emphasized by Hennebelle & Perault (1999 ) the 
radiative time needs to be smaller than the dynamical time for 
the thermal instability to occur. Equivalently, the cooling time 
must be smaller than the typical time of the stirring. 



(y - l)nh 



< L s /v s . 



(19) 



For fixed box size (L s ) and initial temperature, this offers a natu- 
ral explanation why the transition is inefficient at the lowest den- 
sities and the highest amplitudes vs . In our case L s is 20 pc and, 
for typical initial conditions of n = 1 cm" 3 and T = 8000 K, the 
cooling time is t coo \ = 1.7 Myr implying that v s < 11.4km s" 1 . 

5.2. Effect of the compressive modes 



Federrat h et al.| ( |2008| |2010| ) showed on isothermal simulations 
that the ratio of compressible to solenoidal modes in the forced 
velocity field changes the structure of the gas. To understand the 
effect of the nature of the turbulent forcing on the cold gas for- 
mation in the thermal instability frame, we performed 42 simula- 
tions (128 3 pixels) with different values of the projection weight 
f and a fixed initial density no = 1.0 cm" 3 , close to the typi- 
cal WNM density. We voluntarily chose an initial density that 
we showed does not induce CNM formation in the case of the 
natural state of turbulent forcing (f = 0.5, see section 2.2) in 
order to see if a different combination of compressive versus 
solenoidal modes can trigger the phase transition at low density. 
Like previously the velocity field amplitude vs was varied from 
5 to 20 kms" 1 . The velocity field is either purely compressive 
(f = 0), either purely solenoidal = 1) or a mixing of both 
components (0.1 < ( < 0.9). Table [2] summarizes the values 
of Al bs> M± eo , /cnm and cr turb in the stationary regime for the 
42 simulations. The behavior of /cnm and <x tur b with f are also 



presented on Figure 10 



These results show that a turbulent field with a majority of 
solenoidal modes or with an energy equilibrium between com- 
pressible and vortical modes (f > 0.5) cannot trigger the tran- 
sition to CNM with this intermediate value of the initial density 
(Figure 10 left). This behavior is obvious on the representations 
of the integrated density on FigurefTT] the more compressive 



modes the velocity field has (right side), the more cold gas is 
formed. For a WNM initial density of no = 1 cm" 3 , f needs to 
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Table 1. Results of the 128 3 simulations with natural turbulence £=0.5 : the two first columns are the initial density and the 
large scale velocity of the turbulent stirring, the last five are the obtained results: average theoretical and observational Mach 
numbers, average and standard deviation of the cold massive fraction, and turbulent velocity dispersion of the line-of-sight velocity 
component. The average and standard deviation were computed between 20 and 40 Myrs. 
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be smaller than 0.3 to reach 30% of the mass in the cold phase. 
Besides, as previously showed and for the same reason, /cnm 
increases at lower values of the large scale velocity at a given 
value of the spectral weight £ (Fig ITU] left part). 

One can notice that cr turb continually increases with the pro- 
jection weight f (right part of Fig 10), even when no cold gas 
is created in the simulation. The turbulent motions estimated 
with (Tt ur b are thus affected by their repartition in the compres- 
sive and solenoidal modes. <x tur b increases with the energy frac- 
tion injected in the solenoidal modes for the same total amount 
of energy injected. The solenoidal modes are efficient at mix- 



ing and uniformizing the density more efficiently than the com- 
pressive modes, which create higher density contrasts correlated 
with the velocity field (Feder rath et al.|2010 ). The cr tur b compu- 
tation gives more weight to the high density regions and thus to 
the compression regions where the velocity field is close from 
uniform, explaining why cr turb decreases when the compressive 
modes become dominant. 
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Table 2. Results of the 128 3 simulations with initial density n=l cm -3 : the first two columns are the turbulence parameters 
(solenoidal over compressible modes rate and amplitude of the turbulent forcing), the next five are the obtained results: average 
theoretical and observational Mach numbers, average and standard deviation of the cold massive fraction and turbulent velocity 
dispersion of the line-of-sight velocity component. The average and standard deviation were computed between 20 and 40 Myrs. 
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6. High resolution simulations 

All the results presented hereafter are concerning the two 
1024 3 cells simulations after 16 Myr. We selected two different 
sets of initial conditions leading to the formation of 30 to 50% 
of cold gas, a subsonic Mach number and a turbulent veloc- 
ity dispersion between 2 and 3 km s" 1 . The first one (hereafter 
1024N01) has an initial density of 1 cm -3 , a majority of com- 
pressible modes with £ = 0.2, and a large scale velocity vs of 
12.5 km s" 1 . The second one (hereafter 1024N02) has an initial 
density of 2 cm -3 , high enough to allow us to study the natural 
state of the turbulent velocity field with f = 0.5, and a large scale 
velocity vs of 7.5 km s" 1 . 



Both simulations are shown on Figure 12 that represents 
the integrated density along the z-axis. They reveal structures 
in clumps and filaments comparable to observations. Besides, 
1024N02 creates much more structures than 1024N01, suggest- 



Table 3. Values of the theoretical Mach number computed at 
first on all pixels and then only on the pixels with a tempera- 
ture greater than 200 K, of the observational Mach number and 
of the turbulent velocity dispersion for both simulations at high 
resolution. 



Simulation 


Mheo 


M the o(T>200) 


M) bs 


cr turb (kms -1 ) 


1024N01 


0.85 


0.80 


0.47 


2.77 


1024N02 


1.20 


1.20 


0.68 


2.20 



ing that it is more efficient at triggering the transition, which we 
will discuss later. 



11 



Saury et al.: The structure of the thermally bistable and turbulent atomic gas in the local interstellar medium 




0.4 0.6 
spectral weight 



g> 
"c/> 



1.0 




0.4 0.6 
spectral weight 



1.0 



Fig. 10. Behavior of /cnm (on the left) and cr turb {on the right) with the variation of the projection weight f , for the different values 
of the large scale velocity vs (from 5 to 20 km s" 1 ). The dotted lines frame the observational values of /cnm and <T tur b that we want 
to reproduce. 




Fig. 11. Maps of the logarithm of the integrated density along the direction z. The projection weight decreases from the left to the 
right and takes the following values: 0.5, 0.4, 0.2 and 0.0. The initial density and the large scale velocity are fixed to no = 1.0 cm" 3 
and vs = 12.5 km s" 1 . All maps are computed on the same column density range [0.3 10 20 cm -2 , 32 10 20 cm -2 ]. 



6.1. Mach number 

Table [3] summarizes the values of AVeo> computed at first on all 
the pixels, and then only on the pixels with a temperature greater 
than 200 K, At bs and <T tur b. The theoretical Mach numbers show 
that 1024N01 is subsonic while 1024N02 is transsonic. Besides, 
removing the cold pixels of the box (T < 200 K) in the compu- 
tation of Attheo does not change the results, suggesting that the 
turbulent motions are dominated by the WNM motions, as we 
will discuss it while studying the CNM structures. The observa- 
tional Mach number is in both cases clearly subsonic, in agree- 
ment with observations. Lastly, the computed values of <x tur b also 
show that the observed properties of turbulence are well repro- 
duced in both simulations although it is slightly lower than the 
estimated value in a WNM cloud of 40 pc (~ 3 kms" 1 ). This is 
due to the cold structures present in the simulations for which 
the velocity dispersion is much lower. 



6.2. Thermal distribution of the gas 

As expected with the thermal instability at low Mach number, 
the distributions in the P, n, T space present the evidences of a 
multiphasic medium. The temperature histograms (Fig. 13 ) have 
indeed a clear bimodal shape, with a peak in the warm phase 
and one in the cold phase. Similarly the massive distributions of 
the gas in a pressure-density diagram (Fig. 14 ) clearly show two 
preferential zones: the first one at high densities and the second 
one at low densities, suggesting again CNM and WNM. 



However, the properties of the warm gas are slightly different 
in both simulations. If the cold gas peaks at similar temperatures 
in both cases (40 K for 1024N01 and 50 K for 1024N02), it is 
not true for the warm gas: while the WNM of 1024N01 has a 
mean temperature of 7000 K, close from the prediction of the 
two-phases model ( [Field et al.|1969| ), the warm gas in 1024N02 
peaks around 3500 K, two times smaller. This effect can also 
be seen in the pressure-density diagrams (Fig. [14]), where the 
pressure of 1024N02 is smaller than the one of 1024N01 at low 
densities. Therefore, while the gas at low density of 1024N01 
follows the stable branch of the WNM, the gas at low density 
of 1024N02 is shifted below the equilibrium curve where the 
heating is dominant. This pressure difference is confirmed by 
the pressure histograms presented on Figure [l5j where the mean 
pressure of 1024N02 (bottom) is lower than the mean pressure of 
1024N01 (top). On the other hand, the density distributions pre- 
sented on Figure 16 peaks around 0.5 - 0.6 cm -3 in both cases, 
similar to the mean density estimated from observations proper- 
ties. 

Two factors can be considered to explain the presence of 
such a warm gas out of equilibrium. First, raising the initial 
density induces the decrease of the cooling time in the WNM 
(eq.[8]j: t coo \ goes from 1.75 Myr for no = 1.0 cm" 3 to 1.15 Myr 
for no = 2.0 cm" 3 . Therefore, the warm gas at low density lo- 
cated below the equilibrium in 1024N02 is cooled down faster 
than in 1024N01 and does not have time to be heated enough to 
reach the equilibrium curve. Secondly, the thermal equilibrium 
curve considered in the pressure-density diagrams is established 
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for a static gas and does not take into account the energy con- 
tribution of turbulence. Depending on the properties of the tur- 
bulence, WNM gas can stabilized at a temperature significantly 
lower than the prediction of the two-phases model. 

The mass fractions of the different thermal phases have been 
calculated using two different temperature thresholds and one 
density threshold. All results are summarized in Table [4] We 
first note that the amount of CNM created by the simulations 
is 30% for 1024N01 and 50% for 1024N02. These quantities are 
in perfect agreement with the predictions deduced from the 128 3 
simulations and with the observational constraint of Heil es &l 
|Trola nd (2003) who emphasized that 40% of the mass of the H i 
lie in the cold phase. Besides, we estimated in Section [3] the cold 
mass fraction around 44% and the CNM volume filling factor at 
1 % from observations properties of H i. The CNM volume fill- 
ing factor is here also well reproduced with 1% for 1024N01 and 
4% for 1024N02. 

Regarding the warm gas, about 20% of its mass is lying 
between 2000 and 5000 K in both simulations. In 1024N02, 
where the WNM-CNM transition is highly efficient, the unsta- 
ble regime holds the major part of the mass (46% between 200 
and 5000 K) and half of it has a temperature lower than 2000 K; 
in 1024N01, where the WNM-CNM transition is less efficient, 
only 10% are lying in the unstable regime under 2000 K and a 
significant amount of the gas lies in the WNM range beyond 
5000 K (40%). 

If we use a density threshold to separate the phases, we note 
that both simulations hold the same amount of mass of unstable 
gas (~ 30%) and that the difference between the two simulations 
only lies in the CNM and WNM phases : the amount of gas that 
does not transfer to CNM in 1024N01 stays in the WNM phase. 
Besides, in the 1024N01 case, the three thermal phases are 



distri buted like 1/3-1/3-1/3, as observed by Heil es & Tro land 
(2003), while in the 1024N02 case, the mass fraction of WNM is 
lower of a factor 2 and more than 50% of the mass is in the CNM. 

2D density slices with isocontours placed at 200 K (blue) and 
2000 K (red) are presented on Figure [17] The highest tempera- 
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Fig. 13. Temperature histograms weighted by the density: 
1024N01 (solid black line) and 1024N02 (dashed blue line). For 
the 1024N02 simulation, the histogram is normalized by 2 to be 
comparable to 1024N01 



ture threshold has been consciously chosen lower than the usual 
limit at 5000 K in order to see the distribution of the thermally 
unstable gas. The gas lying between 2000 and 5000 K is indeed 
widely distributed and does not show any preferential location. 
On the other hand, the thermally unstable gas below 2000 K lies 
preferentially around the cold structures as sheets around the 
denser regions, in agreement with |GazoT& Kim ( 2 010| ). 



6.3. Pressure range 



Jenkins & T ripp ( 2011) used observations in absorption to study 
the distribution of the thermal pressure in the diffuse, cold neu- 
tral medium. They found that the pressure distribution, weighted 
by the mass density, is well fitted by a lognormal. The pressure 
histograms weighted by the density of our simulations are also 
well represented by lognormal distributions. Figure 15 show the 



3D pressure histograms weighted by the density (black lines), 
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Table 4. Mass fractions for the three thermal phases : cold, warm, unstable range, Mach number and turbulent velocity dispersion 
at 16 Myr for two different criteria in temperature and a density criterion. 



simulation 


/CNM 


/unstable 


/WNM 


CNM volume 


WNM volume 




T<200K 


200K<T<5000K 


T>5000K 


filling factor 


filling factor 


1024N01 


0.30 


0.29 


0.41 


1% 


66% 


1024N02 


0.51 


0.46 


0.03 


4% 


14% 




T<200K 


200K<T<2000K 


T>2000K 






1024N01 


0.30 


0.10 


0.59 


1% 


95% 


1024N02 


0.51 


0.25 


0.24 


4% 


73% 




n > 7.0 cm -3 


0.8 cm -3 < n <7.0cm" 3 


n < 0.8 cm -3 






1024N01 


0.33 


0.32 


0.35 


1% 


75% 


1024N02 


0.54 


0.31 


0.15 


4% 


60% 




1.0 10.0 100.0 1000.0 

density (cc) 




_ 7.00 



10.0 

density (cc) 



1000.0 



Fig. 14. Distributions of mass in a pressure-density diagram 
with on the top the 1024N01 simulation and on the bottom the 
1024N02 simulation. The unit of the 2D histogram is the log- 
arithm of the density (in cm -3 ). The solid blue curve shows 
the thermal equilibrium (cooling equal to heating) with the pro- 
cesses implemented in HERACLES, the dashed-dotted lines are 
the 200 K (green) and 5000 K (yellow) isothermal curves. 
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Fig. 15. Pressure histograms from 1024N01 (top) and 1024N02 
(bottom). The black line is the histogram of the 3 dimensional 
pressure, the blue one is the histogram of the averaged pressure 
along the line of sight. The associated colored dashed lines are 
the lognormal fits. The red dashed-dotted line represent the log- 
normal fit found by Jenkins & Tripp ( 201 1) on observations. The 
lognormal coefficients ( [a\, a{[ - see Eq. 20) are [3.58,0.175] for 
Jenkins & Tripp|([2Ql 1|) For the 3D and averaged pressure, they 



are [3.58, 0.29], [3.58,0.14] for 1024N01 and [3.35, 0.20], [3.38, 
0.09] for 1024N02. 



the averaged pressure weighted by the column density (blue . , . . U1 ^ n . +u , ~ ... 

» i j2 ■% i i j • j_ 'i a' / ■% i j i • n j .1 3-dimensional pressure is not available. Following the definition 

lines), the fitted lognormal distributions (dashed lines) and the ^ Jenkins & Tri (20TTb 

lognormal distribution found by |Jenkins & Tripp ( 2011| ) (red ' »' 

dotted line). The averaged pressure is computed on each line of 

sight and is itself weighted by the density of each cell. By do- & n _ x ^ / _ (logP/k - a\) \ 

ing this, we try to get closer to real observations in which the dlogP/k a ° QX P\ 2 x ^ 
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Fig. 16. Density pdfs for both simulations: 1024N01 (solid black 
line) and 1024N02 (dashed blue line). 



the values found for our lognormal fits are given in the cap- 
tion of Fig. 15 The simulation 1024N01 peaks at the same 
value than their results, meaning 3800 K cm -3 . The pressure 



peak of 1024N02 stays close although slightly lower at 3.35 dex 
(2200 K cm -3 ). These two values lie exactly in t he pressure 
range allowing a two-phase medium according to |Field et aL] 
(1969). As expected, we note that the width of the distributions 
are larger for the 3D pressure histograms than the averaged pres- 
sure histograms in both cases (0.29 and 0.20 compared to 0.14 
and 0.09 respectively). The average along the line of sight indeed 
smoothes the most extreme values. Besides, the distributions of 
1024N01 are larger than 1024N02, which can be explained by 
the lower amplitude of the turbulent stirring used in 1024N02. 
However, the width obtained in these simulations are close to 
the result of Jenkins & Tripp ( 2011| ) and thus are in agreement 
with the observations. 



6.4. Density pdfs 

The density distribution is of great importance in turbulence 
studies. Many simulations of molecular clouds have shown that 
the density distribution of an isothermal gas follows a lognor- 
mal law with a widening that depends on the Mach number 
( Vazquez-Semadeni|| 1994 Passot & Vazquez-Semadeni| |l998; 
|Federrath et al.|2008| ). The density distribution of a two-phases 
gas is however more complex and often presents a bimodal 



shape {jfle nnebelle et aL|2008[| Audit & Hennebelle 



Piont ek & Qstr iker (2005|), |Gazol et al.| ( |2005| ) and |Walcheta 



2005 



2010) 



( 201 1| showed that density pdfs are broadened and loose the bi- 
modal shape when the Mach number increases, as most of the 
gas move to the thermally unstable range. The density distribu- 
tions presented on Figure [16] point the evidence of a transsonic 
biphasic gas with two distinct peaks, one at low densities, 0.4- 
0.5 cm -3 , in agreement with the density estimation in the WNM, 
and the other one around 10-20 cm -3 , in agreement with the den- 
sity in the CNM. 

To compare with studies of molecular clouds, we tried to 
fit lognormals on cold peaks for each simulation. Regarding the 
high density part, it is possible to fit the right tail with a lognor- 
mal, but for both simulations, the lognormals are overly broad. 
The left tail of the lognormals are indeed going to too low den- 
sities to be a good approximation of the CNM density pdf. To 
confirm this, we separated the cells of the box where T < 200K 
and plotted the pdfs of the CNM. They did not show any charac- 
teristic shape, and it was impossible to fit lognormals on them. 



We also did n't find any evidence of a power law tail as Seifried 
et al. (2011 ) did. Therefore, if some lognormal fit is possible on 
the high density part of the pdf, it cannot be used to draw any 
physical conclusion. 

6.5. Power spectra 

Compared with converging flow simulations, the stirring in 
Fourier provides a way to better study the velocity and density 
power spectra and therefore characterize the specific properties 
of the turbulence in a two-phase medium. What is shown here 
is a study of the statistical properties of the density and veloc- 
ity fields in 2D slices along the z-directiorj^] We choose the z- 
component of the velocity and not its norm to stay closer to 
what observations have acces to. For comparison purposes, we 
also ran a 1024 3 cells isothermal simulation, using the same tur- 
bulent stirring than 1024N01 (£ = 0.2 and v s = 12 km s _1 ), an 
initial density of n = 1 cm -3 and a temperature of 8000 K. For 
each simulation, we calculated for each cut along the z-axis the 
density contrast < n > /cr(ri), the maximum of the density and 
computed 2D power spectra of the density and the velocity. We 
present on Figures [18] and [19] the evolution of the contrast, the 
density maximum, and the slopes of the power spectra with z 
and on Figures 20 and 21 the compensated mean power spectra 
k* l?> P(k) of the density and velocity cuts for 1024N01 and for 
the isothermal simulation (these power spectra for 1024N02 be- 
ing similar to those of 1024N01, they are not displayed here). 
The errors bars represent the (lcr) variations on the 1024 cuts. 

As shown by Kim & Ryu ( 2005), the density power spectrum 
of the isothermal and subsonic simulation has a Kolmogorov 
slope of -8/3. We note that the velocity power spectrum also fol- 
lows the Kolmogorov law. The estimation of the slope of the 
power spectrum in the inertial zone of numerically simulated 
turbulence can only be done on a small scale range due to the 
turbulent forcing at large scale and the numerical dissipation at 
small scales. The size of the simulations (1024 3 ) and the use of 
an isothermal simulation are of great help in the estimation of 
the inertial range, since the result is known (P(k) oc /r 8/3 in 2D). 
We deduce from the power spectra of the isothermal simulation 
that the Kolmogorov slope is well reproduced on the decade of 
scales 0.003 < k < 0.03 (or -2.52 < logk < -1.52). All the 
values of slopes fitted on the simulations with thermal instability 
have been computed on this range of scales. 

Regarding the TI simulations, the velocity power spectra 
are slightly flatter from Kolmogorov with slopes of -2.4 for 
1024N01 and -2.3 for 1024N02. On the other hand, the den- 
sity power spectra are far from Kolmogorov with slopes of 
-1.3 for 1024N01 and -1.5 for 1024N02. This is due to the 
much higher density contrast. Its values are indeed around 4 
(< C n0 i >= 4.6 ± 1.4 and < C n0 2 >= 3.8 ± 0.7) while the density 
contrast of the isothermal simulation is almost constant from a 
cut to another and stays close from 0.6 (fig. [18]). 

In supersonic turbulent flows it is generally observed that the 
slope of the density power spectrum flattens with the increase of 
the Mach number. This effect has been reported for isothermal 
(|Kim & Ryu||2005] ) and supersonic, thermally bistable (Gazol 
|& Kim|2 010) flows. This effect is related to the increase of the 
width of the log-normal PDF of the density with Mach num- 
ber (|Vazquez-S emadeni|20 1 2b| ) . It was also noted that while the 
the density power spectrum flattens with the Mach number, the 



2 the choice of the axis is not important as the simulations are 
isotropic due to the turbulent stirring in Fourier space. 
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Fig. 17. Logarithme of a density (cm 3 ) snapshots from the 1024 3 simulations, with contours from the temperature at 200 K (blue 
lines) and 2000 K (red lines). Left: 1024N01, right: 1024N02 



velocity power spectrum stays close from Kolmogorov ( Kritsuk 
|et al.|2007] ), slightly steeper tough. 

We observe a similar behavior on simulations of transsonic, 
thermally bistable turbulence; the density power spectrum is 
much flatter while th e velocity power spectra s tay close to 
Kolmogorov estimate. Hennebe lle & Audit ( |2007| ) emphasized 
that the density fluctuations of the Hi, and thus the density 
contrast, are due to the thermal instability at low Mach num- 
bers while supersonic shocks dominate at high Mach numbers. 
The TI simulations are indeed much more structured at small 
scales than the sub- sonic isothermal one. We also want to point 
out that the power spectrum of the logarithm of the density 
(Pk(log(n)) ), which traces more directly the velocity power 
spectra through the equation of continuity, has a slope very close 
to the Kolmogorov law for both 1024N01 and 1024N02. 

6.6. CNM structures 

Two observational facts often attributed to the self- similarity in- 
duced by the turbulent c ascade in the interstellar medium are the 
Larson||1981[ |Elmegreen & F algarone 



mass-scale (M oc L a 



1996; Roman-Du val et al. 2010) and velocity dis persion-scale 
(cr turb zzo- lpc Z7 - |Larson|1981[|Heyer & Brunt|2004| ) relations ob- 
served in molecular clouds, if CNM clumps are the precursors of 
the molecular matter, it is interesting to see of they follow simi- 
lar laws. To identify CNM clumps in the simulations we selected 
pixels having a density greater than 5 cm -3 . Structures were then 
identified as islands of spatially conn ected pixels. Following the 
approach of Henneb elle et al. ( 2007| ), the size of a structure is 
defined using the largest eigenvalue A\ of its inertial matrix J: 



2 k )dm 



Xjdm 



In = h 



with i, j, k e [1, 2, 3] and i + j ± k 



(21) 

(22) 
(23) 



The Xi,ie 1, 2, 3, are the pixels coordinates and dm the infinites- 
imal mass. One can thus estimate the size of the structure by 
computing L = y/A\/M, where M is its mass. We point out that 



Audit & Hennebelle] ( |2010| ) investigated other choices than the 
largest eigenvalue, as the geometrical mean of the eigenvalues, 
and found very similar results. 

Figure 22 shows the relation between the mass and the size 
of the structures for 1024N01 (this relation for 1024N02 being 
similar to that of 1024N01, it is not displayed here). We found 
slopes of 2.25 ± 0.05 and 2.28 ± 0.03 for 1024N01 and 1024N02 
respectively. This results are in agreement with previous work 
on is othermal simulations ( Krit suk et al.||2007[ |Federrath et al. 
2009 ) and on s imulations including the thermal instability ( Audit 
& Hennebelle 2010). They are also in agreement with the re- 
sults obtained on CO observations by Elmegreen & Falgarone 
(1996 ) who found a value between 2.2 and 2.5 on a large sam ple 
of independent molecular clouds. Recently, Roman-Duval et al. 
( 2010) measured an index of 2.36 on 580 molecular clouds. The 
structure observed here on H i simulations is similar to molecular 
clouds. 

Similarly we investigated the relation between the velocity 
dispersion and the size of the CNM clumps. The velocity disper- 
sion is computed as: 



o-(M) = 



Zp(M-IvqI) 2 
2> 



(24) 



with vo the mean velocity of each structure weighted by the den- 
sity p 



vo 



2> 



(25) 



Figure [23]represents the values of <x(|v|) versus size for the struc- 
tures found in simulation 1024N01 (this relation for 1024N02 
being similar to that of 1024N01, it is not displayed here). 
The fits obtained are 0.70 x L 0Al for 1024N01 and 0.51 x L 0A2 . 



These are compatible with the results of Audit & Hennebelle 
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Fig. 18. Isothermal simulation. Bottom: 2D power spectra slopes 
of density (Pk(n)) and velocity (Pk(v z )) cuts along z. For each 
cut, the plots on the top give the maximum of the density map 
and the density contrast. 
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Fig. 19. Simulations with thermal instability. Bottom: 2D power 
spectra slopes of density (Pk(n)) and velocity (Pk(v z )) cuts along 
z. For each cut, the plots on the top give the maximum of the 
density map and the density contrast. 
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Fig. 21. Mean velocity compensated power spectra Pkk 2 - 67 for 
1024N01 



( |2010| ) who found 0.33 < y < 0.53. They are also in agreement 
with the value s measured on observat ions : 0.37 in Hi clouds 
( |Larson|1979| ) and molecular clouds ( |Larson|1981|), and 0.5 ob- 
tained on a large sample of molecular clouds by Heyer & Brunt| 
(2004). All these values are also close to the estimated value of y 
for subsonic turbulence. On the other hand, the values obtained 
here for cr lp c 0.70 and 0.51 kms" 1 , are lower than the results 
presented by Audit & Henne belle ( |2010| ) who found values be- 
tween 1.2 and 3.3 km s _1 but with a different forcing scheme. The 
turbulence in their simulation is indeed produced by WNM con- 
verging flows with a Mach number close from 1.5. This type of 
forcing is very efficient to create cold structures but induces a 
Mach number in the WNM higher than the one observed in the 
diffuse medium. 

To compare to the work of |Wolfire et al.| ( |2003| ) and to 
the value 0.92 kms -1 estimated in § [5] we computed <x(l) = 
<r(|v|)/Lpc 3 for each CNM structure using the the theoretical 
value y = 1/3 expected for a subsonic turbulence. The his- 
tograms presented on Figure [24] peak around 0.8 and 0.6 kms -1 
for 1024N01 and 1024N02 respectively. These values are in 
agreement with the velocity dispersion at 1 pc we wanted to re- 
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Fig. 22. Structures mass versus their size for 1024N01. 
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Fig. 23. Distribution of the velocity dispersion inside each 
clump versus their size. We overplotted in red the fit of the scat- 
ter 0.70 x L 0A1 



produce in these simulations and with the work of Wolfire et al. 
(2003]). The fact that cr lpc is well reproduced suggests that the 
internal dynamics of the CNM clumps is related to the dynamics 
in the WNM. To confirm this assertion, we computed the mean 
velocity dispersion < <x(|v|) >. We obtained 0.3 and 0.4 km s" 1 , 
much lower than the mean thermal velocities inside the clumps: 
1.7 km s" 1 for 1024N01 and 1.3 km s" 1 for 1024N02. The indi- 
vidual cold structures are therefore subsonic. On the other hand, 
the relative motions between the clumps are greater. We com- 
puted the cloud-cloud velocity dispersion as 



Zo g (VQ- < VP » 2 

Ng 



(26) 



Ng being the number of clumps. We obtained <x cc = 2.3 kms -1 
for 1024N01 and <x cc = 1.9 kms -1 . These values are very close 
from the total velocity dispersions in the boxes (2.8 kms -1 for 
1024N01 and 2.2 km s _1 for 1024N02), suggesting again that the 
clumps motions are related to the motions in the WNM. This 
means that the CNM structures keep the memory of the flow 
in which they were formed and that little dissipation occurs in 
cloud-cloud collisions. We thus conclude that a significant frac- 



Fig.24. Histograms of cr(v z )/L 1/3 for the 1024N01 simulation 
in the solid black line, and 1024N02 in the dashed blue line. 

tion of the observed 21 cm line broadening is due to the relative 
motions of the clumps and is not the result of internal supersonic 
motions. This is in agreement with the work of Heitsch et al. 
( 2005 , 2006 ) who simulated the formation of molecular clouds 
from H i converging flows and reached the same conclusion. 



7. Summary and discussion 

The discussion of the results are divided in two parts. First we 
discuss the physical conditions that favor the formation of CNM 
in the local ISM. Second we describe the results on the physical 
properties of the CNM itself, based on the analysis of the two 
simulations at high resolution (1024 3 ). 



7.1. Formation of the CNM 

1. The first main result of the parameter study on 128 3 sim- 
ulations is that the fiducial conditions of the warm neutral 
medium in = 0.2 cm -3 and T = 8000 K) do not lead to the 
formation of cold gas, whatever the turbulence properties in 
the range we explored. With a slight increase of the den- 
sity (1 cm -3 ), a majority of compressive modes of the stir- 
ring are needed to trigger the transition. A higher increase 
of the initial WNM density (above 1.5 cm -3 ) produces CNM 
very efficiently. A moderate compression effect is therefore 
needed to trigger the transition with either a compressive ve- 
locity turbulent field or with an increase of the WNM density 
(or pressure). To first order the required increase in pressure 
(from Pi to P 2 ) for the WNM-CNM transition could happen 
by the confluence of two flows with a relative velocity A V : 



AV : 



(P2-Pl) 



\pl Pi' 



1/2 



(27) 



Considering n\ - 0.5 cm 3 (typical density in the WNM), 
n 2 = 2.0 cm" 3 (initial density of 1024N02) and T = 8000 K, 
one obtains AV =10.3 km s" 1 . The pressure increase needed 
to place the WNM into suitable conditions to trigger the tran- 
sition can happen with converging flows with a relative ve- 
locity of 10 kms -1 , which is completely feasible in cases of 
supernova? explosions or stellar winds and compatible with 
the idea that turbulence at large scales is produced by con- 
verging flows induced by the star formation activty ([Ferriere 
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2. 



200l||Elmegreen & Scalo|2004||McKee & Qstriker 2007|[de 



Avil lez & Breitschwerdt|2007[|Kim et al.|2010||20iTr 
The second main result is that a sub or trans sonic turbulence 
is required to get CNM gas with properties in agreement with 
observations (/cnm~ 40% and <x tur b~ 3kms _1 ). In oppo- 
sition to an isothermal gas, it is rather inefficient to inject 
supersonic turbulence to produce high density structure in a 
thermally bistable turbulent flow. This result ties up with the 
idea of Vazquez-Semadeni (2012b) and can be explained by 
Eq. [61 when ^yn is lower than t coo \, turbulent compressions 
and oppressions act quicker than the thermal instability and 
prevent its development. On the contrary, when t coo \ < £ dyn , 
the gas has time to cool down while turbulence is compress- 
ing it and to reach the stable branch of the CNM before its 
reexpansion. This argument is in agreement with previous 
studies. Walc h et aT] ( |2011| ) observed that, in the conditions 
of the solar neighborhood, the two phases appear only for 

( |2005| > also gathered that the 



the pressure range allowed by th e two-phases medium ( Field 
et al.|1969 |Wolfire et al.|20 03), even when the initial pres- 



low Mach number. Gazol et al 



ratio ^cooiAdyn increases with the Mach number and that the 
structures become transient and move away from the cold 
branch of the thermal equilibrium. At the same time,|Heitsch| 



et al.| p005| ) deduced from their simulations of H i converg- 



ing flows that the thermal instability is dominant when the 
density is high and the velocity of the flow low. 
3. The third important result is related to the impact of the 
distribution of the turbulent energy between the solenoidal 
and compressive modes. Simulations with moderate initial 
WNM density (no = 1.0 cm" 3 ) and a purely compressive 
forcing field (f = 0) leads to a good fraction of CNM, 
close to what is deduced from observations, but in a very 
long time (> 20Myr). Besides, the Mach number and the 
turbulent velocity dispersion also need more than lOMyr 
to converge to the observed values, suggesting that the 
turbulence takes a long time to be fully developed, namely 
to distribute the energy of the compressive modes towards 
the solenoidal modes that are responsible for the turbulent 
cascade. The convergence and the conversion of the WNM 
in CNM is much faster (~ 1 Myr) in the case of turbulence 
with a natural ratio of solenoidal to compressive modes but 
in this case a slightly higher WNM density (no > 1.5 cm -3 ) 
is required to trigger the transition. Given these timescales, 
we conclude that it is more likely that the formation of CNM 
structures occurs by an increase of the WNM density rather 
than purely compressive turbulent energy injection at lower 
density. 



7.2. Structure of the bistable and turbulent H i 

We summarize here the main results obtained on the H i simula- 
tions at high resolution. 

1. We initially showed that these simulations reproduce well 
the observed physical quantities used as constraints for the 
parameter study : a cold mass fraction around 40% with a 
CNM volume filling factor going from 1 to 4%, a trans sonic 
Mach number and a turbulent velocity dispersion estimated 
around 3 km s" 1 for 40 pc. Furthermore, the mass of the gas 
is distributed in the simulations as it is observed by |Heiles"& 
Troland (2003) between the thermally unstable regime and 



the WNM defined by a density criterion. Finally, the pres- 
sures are in very good agreement with the observations of 
Jenk ins & Tripp| ( |2011| ), for both the mean and the disper- 



sures are higher. 

2. We noted that the thermally unstable gas, defined with tem- 
peratures included between 200 and 5000 K behave in two 
different ways, depending on its temperature. The gas be- 
tween 200 and 2000 K is located at the interface of the cold 
clumps and the warm gas while the gas beyond 2000 K is dy- 
namically connected to the WNM and is widely distributed 
in the box. This wide distribution of t he unstable gas has also 
been noticed by |Gazol & Kim| (2010) who observed that the 
boundaries between CNM and WNM are sharper when the 
Mach number is very low (0.2), due to the unperturbed de- 
velopment of the thermal instability by turbulence. 

3. The density distribution is bimodal, be aring the character- 
istic of a bistable gas, as noticed by |Piontek & Qstriker 



( |2004| ) and |Audit & Hennebelle| ( |2010| ) in a multiphasic 
medium. This characteristic can be seen here because the gas 
is trans sonic. Indeed, a high t urbulence maintains the gas in 
the thermally unstable regime ( [Gazol et al . 2005 ; Walch et al. 



|2011| ) and breaks the bimodality of the distribution that tends 
to be close from a single peak lognormal function. As men- 
tioned by Vazquez-Semadeni (2012a), this is due to the fact 
that for strong turbulence, the gas is not in the thermal equi- 
librium locus of the pressure-density phase diagram but it is 
rather mostly in the thermally unstable range. This might not 
be representative of the Solar neighborhood where the WNM 
is thought to be transsonic, implying a bimodal density PDF. 
In this case we showed that the modelling of the density PDF 
as a lognormal is not convincing, not even the high density 
part where we do not observed a power law tail either. 
4. The velocity power spectra of the simulations with ther- 
mal instability follow the Kolmogorov law (-8/3 in 2D) with 
slopes equals to -2.4 and -2.3 for 1024N01 and 1024N02 re- 
spectively. On the other hand, the density power spectra are 
much flatter (-1.3 and -1.4) and the density contrasts high 
(around 4). These signatures are the direct result of the for- 
mation of dense structure by the thermal instability in a low 
Mach number flow (Hennebelle & Audit 2007). It is inter- 
esting to point out that the power spectra behavior observed 
here is similar to the one seen in supersonic isothermal flows 
where the density power spectrum flattens with the increase 
of the the Mach number ( [Kim & Ryu 2005} while the veloc- 



ity power spectrum stays close to Kolmogorov ( Kritsuk et al. 



2007) or even steeper. 



The mass-scale relation M oc L a found (a = 2.25 and 2.28) 
for clumps with a density high er than 5 cm" 3 are in good 
agreement with observatio ns ([ETmegreen & Falgarone 1996; 
Roman-Duval et al. 20T0|) as well as with previous numer- 
ical studi es, even when different conditions an d algorithm 
are used ([Kritsuk et al.|[2007l |Federrath et ^[20091 |Audit[ 
& Henne belle|2010| ). The y-index of the velocity dispersion- 



sion. It is remarkable that the pressures are stabilized right in 



scale relation (cr oc V) is also well reproduced in both sim- 
ulations with values of 0.41 and 0.42. These values are sim- 
ilar to the ones obtained on numerical simulations ( Audit 
|& He nnebelle 20 10]), to the observ ed in molecular clouds 
( |Larson| 198 1 [|Heyer & Brunt|2004| ), but also to the value ex- 
pected for a subsonic turbulence y = 1/3. These results sug- 
gest that the self- similar density structure observed in molec- 
ular clouds could be inherited from the H i from which they 
formed. 

6. The internal velocity dispersion of the cold structures (0.3- 
0.4 km s" 1 ) indicate that they are clearly subsonic (the ther- 
mal velocity dispersion is 1.3-1.7 km s" 1 ). Like for the 
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WNM, the turbulent velocity dispersion of clumps in- 
creases with scale with a normalisation at lpc of <x(l) = 
0.6 - 0.8 km s" 1 , very close to the value found in the 
WNM. Similarly, the cloud-cloud velocity dispersion is 
1.9-2.3 km s" 1 , close to the total velocity dispersions in 
the boxes (2.2-2.8 km s" 1 ). This implies that the turbulent 
motions inside clumps and the relative velocities between 
them are related to the motions of the WNM from which 
they were formed. Individually, the clouds are subsonic but 
their relative velocities are supe rsonic relative to the CNM 
sound speed (in accordance with Koyama & Inutsuka 2002 ; 
Heit sch et al.||2005|). This result is in comple te agreement 



with the work of [Hennebelle & Audit| ([2007 ) who noticed 
that the CNM velocity dispersion estimated from 21 cm 
spectra is greater than the mean velocity dispersion of the 
CNM structures identified in their simulation. This means 
that a non negligible fraction of the measured velocity dis- 
persion is caused by the relative motions of the clumps along 
the line of sight, suggesting that the observed line broaden- 
ing is likely to be due to the relative clumps motions rather 
than supersonic turbulence. 



8. Conclusion 

We have presented a set of 90 low resolution simulations (128 3 
pixels) in order to study the physical conditions that lead to the 
production of CNM structures out of purely WNM gas, given 
the observational constr aints (velocity dispe rsion, density, pres- 
sure) and the model of Wolfi re et al.| ( |2003| ) of the heating and 
cooling processes at play in the diffuse ISM. These simulations 
show that WNM gas at typical values of density, temperature 
and turbulent motions do not transit to CNM gas. The typical 
volume density of the WNM in the disk lies between 0.2 and 0.5 
cm" 3 . We showed that WNM with initial density under 1 cm -3 
never transit into a cold phase, whatever the turbulent conditions. 
Stronger turbulent motions clearly does not help in producing 
cold gas, in fact we showed that it is the opposite. On the other 
hand a moderate increase of the WNM density (or equivalently 
of its pressure) is very efficient; all simulations with an initial 
density higher than 1.5 cm -3 do transit and the amount of cold 
gas in the simulation increases with the density. When the initial 
density reaches 3 cm -3 , all simulations achieve to create at least 
90% of CNM. These results were obtained with a natural bal- 
ance between compressive and solenoidal modes in the injected 
velocity field. In this case the formation of CNM gas is almost 
instantaneous as soon as the turbulence is fully developed. We 
also explored the case of an intermediate initial density (1 cm -3 ) 
but with varying the ratio of solenoidal to compressive mode in 
the turbulent forcing. We showed that, at this density, the tran- 
sition into CNM only occurs when there is a majority of com- 
pressible modes in the turbulent velocity field but this process is 
much slower (10-20 Myr). 

We conclude that the ISM turbulence cascade on its own 
cannot induce a transition from WNM into CNM at the average 
density of the WNM. The production of CNM is likely to be 
due to turbulent motions of moderate amplitude associated to a 
compressive event of the WNM that could be led by transient 
phenomena such as outflows, supernova explosion or spiral 
density waves. This would trigger the phase transition without 
delay. Once formed, the amount of CNM gas is dynamically 
stable. Therefore there is no need for a constant energy injection 
to maintain the high density contrast of the H i. 



We also have presented a detailed study of the H i based on 
two higher resolution (1024 3 pixels) simulations. To our knowl- 
edge these are the largest simulations to date that include the 
thermal instability and a pseudo-spectral turbulent stirring. The 
bimodal temperature and density histograms show the evidence 
of a bistable medium and the massive fractions of each phase, 
the Mach number of the WNM and the pressure distributions are 
all in agreement with observations. With these simulations we 
confirm a num ber of results of previous studies based on con- 
verging flows ([Heitsch et al.|[2056l |Hennebelle & Audit|[2007l 
|Aud it & Hennebelle 2 010| ) that showed that the structure of the 
CNM and WNM are tightly interwoven: the two media share 
the same velocity fields and the CNM cloud-cloud velocity dis- 
persions is close to the WNM sound speed. While the individ- 
ual CNM structures are clearly subsonic, their relative motions 
are supersonic with respect to their own temperature. In addi- 
tion, thanks to the pseudo- spectral forcing method, the statisti- 
cal properties of the turbulent bi-stable flow could be analysed 
on a larger range of scales than could be done in previous stud- 
ies. We found that for such sub/trans sonic bistable flows, the 
velocity field keeps the properties of subsonic turbulence with 
P(k) oc & -8/3 while the density is highly contrasted, reminis- 
cent of what is observed for supersonic flows and in molecular 
clouds. In the case of the H i simulated here, the joined action of 
the turbulence and the thermal instability allows the formation 
of long-lasting cold and dense structures without need for super- 
sonic motions, providing a favorable terrain for the formation of 
highly contrasted molecular couds. We now intent to use these 
high resolution simulations to create synthetic observations and 
to test and improve some analysis methods used on observations. 

Acknowledgements. Heracle s is available on the following link: 
http://irfu.cea.fr/Projects/Site_heracles/ ( Gonzalez et al. 2007 1. 
The simulations at low resolution have been performed on the cluster Sunnyvale 
at CITA, University of Toronto, and the 1024 3 cells simulations have been 
performed on the cluster Jade, Cines, Montpellier. 
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